APS preprint 



Stochastic Cellular Automata Model for Stock Market Dynamics 

M. Bartolozzi 1 and A. W. Thomas 1 
1 Special Research Centre for the Subatomic Structure of Matter (CSSM), 
University of Adelaide, Adelaide, SA 5005, Australia 
(Dated: February 2, 2008) 

Abstract 

In the present work we introduce a stochastic cellular automata model in order to simulate the 
dynamics of the stock market. A direct percolation method is used to create a hierarchy of clusters 
of active traders on a two dimensional grid. Active traders are characterised by the decision to 
buy, &i(t) = +1, or sell, <Ji(t) = —1, a stock at a certain discrete time step. The remaining cells are 
inactive, <7j(i) = 0. The trading dynamics is then determined by the stochastic interaction between 
traders belonging to the same cluster. Extreme, intermittent events, like crashes or bubbles, are 
triggered by a phase transition in the state of the bigger clusters present on the grid, where almost 
all the active traders come to share the same spin orientation. Most of the stylized aspects of the 
financial market time series, including multifractal proprieties, are reproduced by the model. A 
direct comparison is made with the daily closures of the S&P500 index. 
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I. INTRODUCTION 

Since the successful application of the Black-Scholes theory for option pricing in 1973 
more and more physicists have been attracted by the idea of understanding the behaviour of 
the market dynamics in terms of complex system theory, where self-organized criticality 
and stochastic processes play important roles. The aim of the microscopic models 

proposed so far (for general reviews 0,Q]) is to reproduce some stylized facts concerning 
the temporal fluctuations of the price indices, P(t). In particular, the logarithmic price 
returns 

R(t)=lnP(t + l)-lnP(t), (1) 
and the volatility, defined in the present work as 

v(t) = \R(t)l (2) 



have been studied extensively |J| from an empirical point of view. The results have shown 
that while long time correlations are present in the volatility, a phenomenon known as 
volatility clustering, they cannot be found in the time series of returns. Moreover the latter 
show an intermittent behaviour that recalls in some aspects hydrodynamic turbulence 0, 
llfll . lri| . characterized by power law tails in the probability distribution function (pdf). 
Microsimulations have demonstrated that this kind of behaviour can originate both as a 



stochastic process with multiplicative noise |l2, 



16 



13 . I4I ] and as a percolation phenomenon 
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In order to reproduce these features of real markets we introduce a stochastic cellular 

automata model, representing an open market. That is, a market where the number of 

active traders, defined as cells with spin state different from 0, namely (Xj(t) = ±1, evolve 

in time according to a percolation dynamics. The percolation dynamics is chosen in order 

I I 

to simulate the herding behaviour typical of investors [28|. According to this, active traders 
gather in clusters or networks where, following a stochastic exchange of information, they 
formulate the trading strategy for the next time step. The results obtained by the simulations 



are then compared with the time series of daily closures of the S&P500 index |4 , 



28, 



29, 



30 



over a period of about 50 years 



Moreover, recently, the fractal properties [18] of the price fluctuations have also been 
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investigated for different markets 
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22 1 . A common feature found in these studies 



is the existence of a nonlinear, multifractal spectrum that excludes the possibility of efficient 
market behaviour j^j. The origin of the muftifractality of in the financial time series has 



also been at the center of discussions 
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271 ] . In this paper we consider the multi- 



fractal spectrum of the price fluctuations as a stylized fact of the market time series without 
addressing any question about the underlying process able to generate it. The multifractal 
spectrum is used as a further test for our model. 

A parallel between multifractal and thermodynamical formalism has also been investi- 
gated. We found, in agreement with the previous work of Canessa [45], that the analogue 
specific heat can provide a good tool to characterize intermittency, that is financial crashes 
or bubbles, from a thermodynamics-equivalent point of view. 

II. THE MODEL 



In the present work we simulate the financial market dynamics via a stochastic cellular 
automata model. The agents of the market are represented by cells on a two dimensional 
grid, 512x128. The ith agent at the discrete time step t is characterized by three possible 
states or spin orientations, o~i(t) = 0, ±1. The value o~i(t) = +1 is associated with the 
purchase of a stock while (7j(t) = — 1 with selling. The former states are called active. The 
cells with spin value crj(i) = are inactive traders. The active traders herd in networks or 
clusters via a direct percolation method related to a forest fire model j3] . The information 
carried by the active traders, that is their spin state, is shared with the other members 
of the cluster. The percolation dynamics allow a time dependent herding behaviour and 
the market can be interpreted as an open system not bounded by conservation laws. The 
clustering process will be discussed in detail in the next subsection. 

The trading dynamics is instead related to the synchronous update of the spins of the 
active traders, ruled by a stochastic exchange of information between them, similar to a 
random Ising model 3, Q]- A particular feature of the present simulai 
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3, 



but 



information is not spread all over the grid, as in other multi-agent scenarios 
it is limited by the clusters of interaction previously defined. The mechanism for the spin 
dynamics is explained in "stochastic trading dynamics" subsection. 
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A. Percolation Clustering 



One of the aims of our cellular automata model is to reproduce the herding behaviour 
of active traders We refer to herding behaviour as the tendency of people involved 

in the market to aggregate in networks or clusters of influence. The traders then use the 
information obtained by their network in order to formulate a market strategy. Even if 
the topological structure of these networks of information is not important, since several 
kinds of long range interaction are available nowadays ^J, the number of connections for 
each trader must be, in any case, finite and not extended over the whole market. In this 
framework a direct percolation method is used to simulate herding dynamics between active 
traders. If we assume that the neighbours of influence are those of von Neumann (up, down 
, left, right), the percolation is fixed by the following parameters: 

Ph- the probability that an active trader can turn one of his inactive neighbours into an 
active one at the next time step, (Tj(t) = — > o~i(t + 1) = ±1. This simulates the fact that 
certain information possessed by a trader may induce a potential trader to join the market 
dynamics. 

Pd- the probability that an active trader diffuses and so becomes inactive, o~i{t) — ±1 — > 
o~i(t + 1) = 0, due to each of his inactive neighbour. This mimics the fact that only traders 
at the borders of a network, that is the weaker links, can quit the market. 

p e : the probability that a non trading cell spontaneously decides to enter the market 
dynamics, o~i(t) = — > ai{t + 1) = ±1. 

The values of the adimensional parameters, pt, Pd and p e , influence the stability of the 
system and the percentage of active traders on the grid. In order to test different market 
activities we fix the values pd = 0.05 and p e = 0.0001 while we tune the parameter ph.. At 
the beginning of the simulation the grid is loaded randomly with a small percentage of active 
traders and then the system is permitted to evolve according to the previous rules. If we 
are in a stable range of the parameter ph, after a transient period that depends both on the 
parameter values and the initial number of active cells, the number of active traders on the 
grid begins to fluctuate around a certain average, as shown in Fig. ^ (Top). The market 
can be considered open since the number of agents changes dynamically in time. In this 
regime, the competition between herding and diffusion produces a power law distribution of 
the cluster size, as shown in Fig. E] (Bottom), p(S) ~ S~ x , where S is the cluster dimension, 
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FIG. 1: Top: Different values of the parameter H produce different activities of traders on the 
grid. Bottom: Cluster size distribution for ph = 0.0485 at t = 9000. 



defined as the number of active cells belonging to the same cluster, and A > 0, creating a 
hierarchy of networks. This hierarchy is necessary if were are to take into account a real 
aspect of the market, namely that different traders also have different trading powers. A 
reasonable assumption is that people having a larger number of sources of information, so 
belonging to greater clusters, can be associated with professional investors that, most likely, 
are able to move a greater amount of stocks compared to the occasional investor. Using this 
assumption we are able to define a proper weight for the trading power of different cells, as 
we will discuss in the next subsection. 

A similar percolation model has also been used to reproduce some statistical and geo- 
metrical features of solar activity 0, 2^. 



B. Stochastic Trading Dynamics 

The dynamics of the spins of the active traders, of (i) = ±1 for i = 1, N k (t) (where the 
superscript k, from now on, refers to the kth cluster of the grid configuration at time step 
t) follows a stochastic process that mimics the human uncertainty in decision making [13]. 
Their values are updated synchronously according to a local probabilistic rule: cr 4 fc (t+l) = +1 
with probability and af(t + 1) = —1 with probability 1 — p*. The probability p\ is 
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determined, by analogy to heat bath dynamics with formal temperature fc&T = 1, by 



where the local field, /■ (t), is 



^) = ^E44W + ^ ( 4 ) 



The A^-(t) are time dependent interaction strengths between agents and h k (t) is an external 
field reflecting the effect of the environment The interaction strengths and the external 
field change randomly in time according to A k j(t) = A^ k {t) + arjij(t) and h k (t) = hQ{t). 
The variables £ k (t),r]ij(t),(i{t) are random variables uniformly distributed in the interval 
(-1,1) with no correlation in time or space. The measure of the strengths of the previous 
terms, A, a and h, are constant and common for all the grid. 

In this contest the dynamics of the price index, P(t), can be easily derived if we assume 
that the index variation is proportional to the difference between demand and supply, 

dP 

- oc X P, (5) 
and using a weighted average for the orientation of the spins, 

N cl (t) N k (t) 

= £iv fc (tK fc (t), (6) 
fc=i i=i 

where N c i(t) and N k (t) are, respectively, the number of clusters on the grid and the size 
of the kth cluster, while (3 is a normalization constant. The relation © follows from the 
assumption of proportionality between the financial power of an active cell and the size of 
the cluster to which it belongs. A justification of this relation is that professional investors, 
able to move larger amounts of capital, are more likely to be linked with a large number 
of other investors than the occasional traders. From (jSJ) we find that the logarithm of the 
price returns that is the fundamental quantity we aim to model, is proportional to the 
average orientation defined in (0), R{t) oc x{t). 



III. NUMERICAL RESULTS 



We compare the results of our model with the Standard & Poor 500 (S&P500) index 
that is one of the most widely used benchmarks for US equity performance. The database 
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FIG. 2: Top: Normalized logarithmic returns for the S&P500. Bottom: Time series of returns 
reproduced with the simulation with p^ = 0.0493, A = 1.8 and h = 0. 

analyzed is composed by of the daily indices from 3/1/1950 to 18/7/2003 for a total of 13468 
data. The time series of the index prices, P(t), is converted into the logarithmic returns ([TJ 
and then normalized over the time interval, T, 



r{t) 



R{t) - (R(t))r 



(7) 



a(R(t)) 

where (• • -)t denotes the temporal average and a is the standard deviation. In this way all 
the sets of returns will have zero average and a = 1 . 

Before discussing the results we briefly describe the behaviour of the cellular automata 
with respect to changes of the parameters. Regarding the percolation, the herding parameter 
Phi as previously seen, determines the concentration of active traders on the grid. For a low 
concentration the active traders will be distributed in small clusters and the information 
will be extremely split up on the grid. In this case large coherent events will be more rare. 
A higher concentration, where clusters of the order of a thousand cells are present, allows 
large events to occur with higher frequency. In fact, the phase transition of large clusters 
can easily trigger a crash or a bubble in the market. Time series related to a higher herding 
parameter are therefore more intermittent. We can then infer that a market made by small 
groups of traders behaves like a noisy market, while big crashes or bubbles must necessarily 
be related to the interaction of large networks and so to a kind of crowd behaviour. For 
most of the simulations we fixed the parameter ph = 0.0493 in order to have bigger clusters, 
with S max {t) of the order of 2500 cells. 
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The strength A also plays an important role in the trading dynamics. With ph fixed, 
this parameter is related to the intermittency of the system. Both for large values of the 
activity (A > 10) and for A — > we observe an approach of the pdf toward a Gaussian-like 
shape. That is very, large fluctuations become more and more rare, and A can be regarded 
as a temporal scale for the system, similar to the activity parameter in the Cont-Bouchaud 
model 16j. In spite of this some large fluctuations can be still identified. This is probably 
one of the main differences between the Cont-Bouchaud model and the present. In fact 
Monte Carlo simulations of the former show that an increase of the activity brings a 
rapid convergence toward a Gaussian distribution because a large number of clusters are 
trading at the same time following a random procedure of decision making jl^l . 17 1: there 
are no clusters that can influence the market more than others and so the resulting global 
interaction is noise-like. In our model fluctuations are always allowed because of the heat 
bath dynamics. Clusters of active traders can always be subjected to phase transitions, 
independently of the state of other clusters, creating a displacement between demand and 
supply. 

In order to reproduce the behaviour of financial time series we work with 1.5 < A < 2.5. 
In this range there is a nonlinear dependence of the intermittency on A. The value of the 
parameter a is fixed by the relation a = 2A, already used in the work of Krawiecki et 
al. However, we note that the ratio a/ A is actually not essential in reproducing the 

intermittency found in real data. Even the parameter h does not have a central role in 
the simulation, as long as h « 1. In fact, the percolation dynamics of the active traders 
already introduces a natural noise. In most of the runs we simply set h = 0. 

Now we discuss the results of the cellular automata. In Fig. |21 (Top) and Fig. (Bottom) 
the normalized logarithmic returns of the S&P500 and of the simulation are shown, respec- 
tively. The average number of active traders, in the stable regime, is ~ 16000, as shown in 
Fig. IU(Top). The model reproduces the intermittent behaviour of the S&P500 time series, 
as expressed by the leptokurtic pdf of FigH3 The tails of the distribution follow a power law 
decay, reflecting the fact that large coherent events, far from the average, are likely to occur 
with a frequency higher than expected for a random process (where the shape would be a 
Gaussian). These large events are related to financial crashes or bubbles of the market and, 
in our model, to a phase transition in the spin state of large networks of active traders, as 
we will discuss further on. From a power law fit, p(r) ~ r~ 1-7 (for \r\ > 2), we find 7 ~ 3 
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FIG. 3: Probability distribution function for the S&P500 computed with daily data from 31/1/1950 
to 18/7/2003 and the model. A Gaussian is also plotted for comparison. The parameters used are 
p h = 0.0493, A = 1.8 and h = 0. 

for both the S&P500 and the model, confirming the good agreement between the two. 

The problem of finding the best distribution describing the price returns is a very impor- 
tant issue from a practical point of view jj]. The standard Black-Scholes theory for option 
pricing Q, U, 0] assumes that the returns are normally distributed. This fact has been proven 
to be empirically false, as shown also in Fig. El (see Ref. j4j for a general reference). Finding a 
more appropriate distribution would be an important improvement in this field of research*. 

In order to understand the trading dynamics of the automata we also show two snapshots 
of the grid configuration, the first during a normal session, FiglU (Top), and the second 
during a crash, Fig0] (Bottom). During the normal session the orientations of the spins 
are distributed uniformly over the various clusters and there is no sharp difference between 
demand and supply. The situation is different during a crash. In this case the clusters at the 
top of the hierarchy, the bigger ones, play a fundamental role. In fact they undergo a phase 
transition where the greatest part of their spins share the same orientation. The capacity of 
the clusters to generate a coherent orientation of the spins, and hence of their trading state, 
can be interpreted in terms of a multiplicative noise process , where the collective 

synchronization arises as a result of the randomly varying interaction strengths between 



* An intriguing framework of investi gati on has been provided by the non- extensive statistical mechanics 
proposed by Tsallis EE EE EE EE Elf • A more complete discussion on this important topic is beyond 
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the scope of this paper and will be discussed in future work. 
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FIG. 4: Top: snapshot of the grid during a normal trading period. The black cells are the buyers 
while the orange ones are the sellers. The parameters used in this simulation are ph = 0.0493, 
A = 1.8 and h = 0. Bottom: the same simulation during a crash. Large clusters of sellers are 
indicated by arrows. 

agents. The peculiarity of our model is that crashes or bubbles (sudden price changes) are 
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but rather to a phase transition 



related not to a phase transition of the whole market 
in one or more of the larger clusters that have a greater influence on the trading session. 
This behaviour is probably closer to the real market where the synchronization of trading 
opinion is more likely to happen between large groups of traders than over the whole market. 

The temporal correlations of the logarithmic returns and of the volatility are investigated 
via the autocorrelation function, defined as 

T-T 



c {r) = £>(t + r)z(t), 



(8) 



t=i 



where Tis the length of the time series and r is a time delay for the normalized variable x(t). 
The results for both the model and the S&P500 are shown in Fig|3](Top) and FigEl (Bottom), 
respectively. While the temporal correlation for the returns is lost almost immediately, the 
volatility manifests a slow decay in time, related to the phenomenon of volatility clustering. 
The previous temporal dependencies have been found in both real data and in the simulation. 
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FIG. 5: Top: Autocorrelation function for the price returns. Bottom: Autocorrelation function 
for the volatilities. In both the graphs the parameters used for the model are ph = 0.0493, A = 1.8 
and h = 0. 

IV. MULTIFRACTAL ANALYSIS 



It is also worth pointing out that financial time series present an inherent multifractal 
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221 has been addressed 



ity [lsj]. In the past few years the work of many authors (19, 
to the characterization of the multifractal properties of financial time series, and nowadays 
multifractality can be considered as a stylized fact. In order to study the multifractal prop- 
erties of our model we use the generalized Hurst exponent [35], H(q), derived via the q— order 
structure function, 

S q (T) = (\x(t + T)-x(t)\*) T OCT* B <!'\ (9) 

where x(t) is a stochastic variable over a time interval T and r the time delay. The gener- 
alized Hurst exponent, defined in @, is an extension of the Hurst exponent, H, introduced 
in the context of reservoir control on the Nile river dam project, around 1907 jl8L l36j| . This 
technique provides a sensitive method for revealing long-term correlations in random pro- 
cesses. If H(q) = H for every q the process is said to be monofractal and H is equivalent to 
the original definition of the Hurst exponent. This is the case of simple Brownian motion 
or fractional Brownian motion. 

If the spectrum of H(q) is not constant with q the process is said to be multifractal. From 
the definition © it is easy to see that the function H(l) is related to the scaling properties 
of the volatility. By analogy with the classical Hurst analysis, a phenomenon is said to be 
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FIG. 6: Multifractal spectra for the S&P500 in the period 31/1/1950 to 18/7/2003 (top) and the 
model (bottom). For the latter p h = 0.0493, A = 1.8 and h = 0. 

persistent if H(l) > 1/2 and antipersistent if H(l) < 1/2. For uncorrelated increments, as 
in Brownian motion, H(l) = 1/2. In Fig. [HI a comparison is shown between the multifractal 
spectra of the model and the S&P500 obtained from the prices time series. It is clear that 
both processes have a multifractal structure and the price fluctuations cannot be associated 
with a simple random walk as in the classical efficient market hypothesis 3| . 

The multifractality of the time series can also be discussed in terms of thermodynamic 
equivalents, according to multifractal physics 2, Q- 44, Q]- m this approach we divide 
the time series, x(t) , for t = 1,...,L into TV equal sub-intervals. Then we can write the 
following measure for each of these, 



\x(t + r) - x(t)\ 
Y^ =1 \x(t + r)-x(t)\ 



(10) 



with i — 1, M and r the time delay. The quantity yUj(r) can be viewed as a normalized 
probability measure. The corresponding generating function is given by 



Xq 



i=l 



that is an analogue of the partition function in thermodynamics. According to Ref. |4 
scaling exponent \ q is directly related to the generalized multifractal dimension D q jl8] 



(11) 
the 



Xq = fa ~ l ) D q- 



(12) 
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If we consider Xq as the free energy of our system, Eq.(|12|) provides a link between 
the classical thermodynamical formalism and multifractality. Assuming q as an equivalent 



temperature, we can define an analogue specific heat 
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46|, 



a, 



w (13) 

The previous equation for C q can also be written in terms of singular measure formal- 
ism [I?}]. In this framework we define a measure, e(l,t), as 



e(l,f) 



|ar(t + 1) - x(t)\ 



(14) 



We can then generate a series of measures on shorter intervals of length 5, e(S, I) where 
5 is an integer power of two and I is the index of the first point of the sub segments at that 
resolution. The average measure in the interval [I, I + 8] is 



1 



l+S-l 



e(S,l) = - £ 6(1, r 

6 l*=l 



(15) 



for I — 0, L — 5. In this case we have a scaling property for the ensemble average with 
respect to the scale 5: 



{e(6,l)) oc5 



-K a 



(16) 



In a multifractal process the exponent K q is a nonlinear function of orelated to the inter 



mittency of the time series and to the generalized dimension via 



{q-l)D q = q-l-K q . 
From (fT3|) and (fTTj) we have that [46], 



(17) 



a, 



d 2 K n 



(18) 



q dq 2 ' 

Following Ref. [45j we have found the analogue specific heat for both the S&P500 and for the 
price time series generated by our model, see Fig|7| For r = 1 we observe a double-humped 
shape for both the model and the empirical data. If we take longer time delays the shoulder 
on the right hand side disappears, leaving only a sharp peak, similar to a first order phase 
transition around q = —1.5. The results are in agreement with the analysis of Canessa 
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FIG. 7: Top: analogue specific heat for the S&P500 from 31/1/1950 to 18/7/2003 for two different 
time delays, namely r = 1 and r = 100. A sharp peak is clearly visible around q = —1.5. The 



second peak on the right hand side disappears increasing the tempora 
been computed for the logarithm of the price using the algorithm in 



delay. The C q curves have 



47j for K q . Bottom: analogue 



specific heat for the model with parameters = 0.0493, A = 1.8 and h = 0. The double humped 
shaped for small temporal delays is visible also here. 

and recall the Hubbard model for small to intermediate values of the local interaction [43:]. 
As also suggested in j^J, the second peak is due to the large fluctuations at small scales, 
that is crashes and bubbles. Increasing the time delay means that the fluctuations tend to 
be smoothed and the time series of returns approach a noise-like regime. For this reason 
the analogue specific heat shapes for r = 100 are basically indistinguishable. From this 
argument we can interpret the analogue specific heat, and in particular the second peak for 
low time delays, as a way to characterize crashes, or in general the degree of intermittency 
in a time series. The difference in the shapes and heights of the shorter peak for r = 1 is 
due to a slightly different correlation of the fluctuations in the two time series. Moreover, 
in the model we can link the shorter peak to the physical phase transition in the spin state 
of a network of traders. 
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V. CONCLUSIONS 



In this paper we have introduced a stochastic cellular automata model for the dynamics 
of the financial markets. The main difference between our model and other stochastic simu- 
lations based on spin orientation of agents 12 , k| is the temporal evolution of the networks 



of interaction and therefore the concept of an "open" market. The active traders follow a 
direct percolation dynamics in order to aggregate in networks of information. This makes 
our simulation, even if still a raw approximation, surely closer to the real market, where no 
conservation rules for the number of agents can be claimed. Crashes and bubbles can be 
interpreted as a synchronization of the spin orientation of the more influential networks in 
the market. Moreover, the introduction of a limitation in the number of interacting agents 
reduces drastically the number of computations on the grid per time step. In a system where 
all the agents interact with each other this number goes like N%, being N a the number of 
active agents, while in our model, considering the distribution of the clusters, it is easy 
to see that it goes like N%~ x . The value of A found for several herding parameters, ph, is 
A ~ 0.6 1.4, so that the computational cost is much lower for the present model. This gives 
one the possibility to simulate the market using a very large range of agents. The model 
is able to reproduce most of the stylized aspects of the financial time series, supporting 
the idea that crashes and bubbles are related to a collective synchronization in the trad- 
ing behaviour of large networks of traders, where the information is exchanged according 
stochastic interaction between them. 
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